tester_2.0
William Norfolk
8/27/2020
Load in Libs
library(broom)## Warning: package 'broom' was built under R version 3.6.3library(ggpubr)## Warning: package 'ggpubr' was built under R version 3.6.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.6.3library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.6.3## -- Attaching packages ------------------------------------------------------------------- tidyverse 1.3.0 --## v tibble 3.0.3 v dplyr 1.0.0
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## v purrr 0.3.4## Warning: package 'tibble' was built under R version 3.6.3## Warning: package 'readr' was built under R version 3.6.3## Warning: package 'dplyr' was built under R version 3.6.3## Warning: package 'forcats' was built under R version 3.6.3## -- Conflicts ---------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(mgcv)## Warning: package 'mgcv' was built under R version 3.6.3## Loading required package: nlme## Warning: package 'nlme' was built under R version 3.6.3##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.library(plotly)## Warning: package 'plotly' was built under R version 3.6.3##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutRead in raw data from RDS.
raw_data <- readRDS("./n1_n2_cleaned_cases.rds")
final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
rename(Facility = wrf) %>%
mutate(Facility = recode(Facility,
"NO" = "WRF A",
"MI" = "WRF B",
"CC" = "WRF C"))Data split for main figure
#split data to layer for main plotly figure
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>%
select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke)) %>%
group_by(date) %>% summarise_if(is.numeric, mean)
#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')x_top <- max(only_background$new_cases_clarke, na.rm = TRUE)
#creates the two panels of the main figure
p1 <- only_background %>%
plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~new_cases_clarke,
type = "bar",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Daily Cases: ', new_cases_clarke),
alpha = 0.5,
name = "Daily Reported Cases",
color = background_color,
colors = background_color,
showlegend = FALSE) %>%
layout(yaxis = list(title = "Athens Daily Cases", showline=TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#renders the main plot layer two as seven day moving average
p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke,
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
name = "Seven Day Moving Average Athens",
line = list(color = seven_day_ave_color),
showlegend = FALSE)
#renders the main plot layer three as positive target hits
p2 <- plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n1,
symbol = ~Facility,
marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n2,
symbol = ~Facility,
marker = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(yaxis = list(title = "SARS CoV-2 Copies/L",
range = c(3, 8), showline = TRUE,
type = "log",
automargin = TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#adds LOD line
p2 <- p2 %>% plotly::add_segments(x = as.Date("2020-03-14"),
xend = ~max(date + 10),
y = 3571.429, yend = 3571.429,
opacity = 0.35,
line = list(color = "black", dash = "dash")) %>%
layout(annotations = list(x = as.Date("2020-03-28"), y = 3.8, xref = "x", yref = "y",
text = "Limit of Detection", showarrow = FALSE))
p1## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.## Warning: Ignoring 1 observations p2## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.#create smoothing data frames
#n1
smooth_n1 <- only_n1 %>% select(-c(Facility)) %>%
group_by(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke) %>%
summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
ungroup() %>%
mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
mutate(target = "N1")## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X7_day_ave_clarke' (override with `.groups` argument)#n2
smooth_n2 <- only_n2 %>% select(-c(Facility)) %>%
group_by(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke, cases_per_100000_clarke) %>%
summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
ungroup() %>%
mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
mutate(target = "N2")## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X7_day_ave_clarke' (override with `.groups` argument)#loess at 0.80 span (same as ggplot default)
sumfit_1 <- loess(log_sum_copies_L ~ new_cases_clarke, data = smooth_n1, span = 0.8)
sumfit_2 <- loess(log_sum_copies_L ~ new_cases_clarke, data = smooth_n2, span = 0.8)
#generalized additive model
sumfit_3 <- gam(smooth_n1$log_sum_copies_L ~ smooth_n1$new_cases_clarke)
sumfit_4 <- gam(smooth_n2$log_sum_copies_L ~ smooth_n2$new_cases_clarke)
#generalized linear model
sumfit_5 <- glm(log_sum_copies_L ~ new_cases_clarke, data = smooth_n1)
sumfit_6 <- glm(log_sum_copies_L ~ new_cases_clarke, data = smooth_n2)#make a dataframe to average both genes, remove target column, group by date, then average
both_clean <- full_join(smooth_n1, smooth_n2) %>%
select(-c(target)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup()## Joining, by = c("date", "cases_cum_clarke", "new_cases_clarke", "X7_day_ave_clarke", "cases_per_100000_clarke", "sum_copy_num_L", "log_sum_copies_L", "target")combofit_1 <- loess(log_sum_copies_L ~ new_cases_clarke, data = both_clean, span = 0.8)
combofit_2 <- lm(log_sum_copies_L ~ new_cases_clarke, data = both_clean)p4 <- plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~log_sum_copies_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Copies/L: ', round(log_sum_copies_L, digits = 2)),
data = both_clean,
marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_lines(x = ~date, y = predict(combofit_1),
data = both_clean,
hoverinfo = "text",
text = NULL,
showlegend = FALSE,
line = list(color = '#1B9E77'))%>%
layout(yaxis = list(title = "Log SARS CoV-2 Copies/L",
showline = TRUE,
automargin = TRUE))
p4#make a plot that stacks the averge of both genes
test1234 <-
plotly::subplot(p4,p1, # plots to combine, top to bottom
nrows = 2,
heights = c(.6,.4), # relative heights of the two plots
shareX = TRUE, # plots will share an X axis
titleY = TRUE
) %>%
# create a vertical "spike line" to compare data across 2 plots
plotly::layout(
xaxis = list(
spikethickness = 1,
spikedash = "dot",
spikecolor = "black",
spikemode = "across+marker",
spikesnap = "cursor"
),
yaxis = list(spikethickness = 0)
)## Warning: Ignoring 1 observationstest1234#add trendlines
#extract data from geom_smooth
#n1 extract
#*********************always update the n = TOTAL NUMBER OF DAYS**************************
cc <- ggplot(smooth_n1, aes(x = date, y = log_sum_copies_L)) +
stat_smooth(aes(outfit=fit_n1<<-..y..), method = "loess", color = '#1B9E77',
span = 0.8, n = 85)## Warning: Ignoring unknown aesthetics: outfit#n2 extract
oo <- ggplot(smooth_n2, aes(x = date, y = log_sum_copies_L)) +
stat_smooth(aes(outfit=fit_n2<<-..y..), method = "loess", color = '#1B9E77',
span = 0.8, n = 85)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#n1
cc## `geom_smooth()` using formula 'y ~ x'
fit_n1## [1] 10.78704 10.99953 11.20416 11.40084 11.58946 11.76994 11.94217 12.10607
## [9] 12.26172 12.40933 12.54893 12.68057 12.80427 12.92009 13.02806 13.12821
## [17] 13.22060 13.30524 13.38219 13.45149 13.51317 13.56726 13.60743 13.62868
## [25] 13.63314 13.62292 13.60013 13.56689 13.52529 13.47747 13.42552 13.37157
## [33] 13.31772 13.26608 13.21878 13.17791 13.12492 13.04673 12.95372 12.85623
## [41] 12.76463 12.68928 12.64055 12.60595 12.56704 12.52519 12.48179 12.43823
## [49] 12.39589 12.35615 12.32041 12.29004 12.26643 12.25096 12.24503 12.25001
## [57] 12.26728 12.30137 12.35238 12.41512 12.48441 12.55504 12.62183 12.67958
## [65] 12.73378 12.79317 12.85739 12.92611 12.99898 13.07566 13.15581 13.23908
## [73] 13.32514 13.41364 13.50424 13.59659 13.69036 13.78520 13.88104 13.97827
## [81] 14.07722 14.17822 14.28161 14.38772 14.49690#n2
oo## `geom_smooth()` using formula 'y ~ x'
fit_n2## [1] 10.64237 10.87620 11.10220 11.32036 11.53061 11.73292 11.92725 12.11355
## [9] 12.29190 12.46242 12.62510 12.77994 12.92694 13.06609 13.19741 13.32087
## [17] 13.43649 13.54426 13.64418 13.73625 13.82046 13.89682 13.96102 14.00972
## [25] 14.04433 14.06628 14.07701 14.07793 14.07047 14.05606 14.03612 14.01207
## [33] 13.98535 13.95738 13.92958 13.90339 13.86050 13.78894 13.69959 13.60336
## [41] 13.51112 13.43377 13.38219 13.34281 13.29590 13.24320 13.18640 13.12722
## [49] 13.06739 13.00861 12.95259 12.90106 12.85572 12.81829 12.79049 12.77402
## [57] 12.77061 12.77976 12.79847 12.82469 12.85635 12.89142 12.92783 12.96354
## [65] 13.00070 13.04279 13.08964 13.14111 13.19706 13.25736 13.32184 13.39037
## [73] 13.46281 13.53901 13.61883 13.70212 13.78874 13.87854 13.97153 14.06787
## [81] 14.16772 14.27118 14.37840 14.48952 14.60465#assign fits to a vector
n1_trend <- fit_n1
n2_trend <- fit_n2
#extract y min and max for each
bb<- ggplot_build(cc)$data## `geom_smooth()` using formula 'y ~ x'bb<- as.data.frame(bb)
n1_ymin <- bb$ymin
n1_ymax <- bb$ymax
qq <- ggplot_build(oo)$data## `geom_smooth()` using formula 'y ~ x'qq <- as.data.frame(qq)
n2_ymin <- qq$ymin
n2_ymax <- qq$ymax
#reassign dataframes (just to be safe)
work_n1 <- smooth_n1
work_n2 <- smooth_n2
#fill in missing dates to smooth fits
work_n1 <- work_n1 %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_n1 <- work_n1$date
work_n2 <- work_n2 %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_n2 <- work_n2$date#create a new smooth dataframe to layer
smooth_frame_n1 <- data.frame(date_vec_n1, n1_trend, n1_ymin, n1_ymax)
smooth_frame_n2 <- data.frame(date_vec_n2, n2_trend, n2_ymin, n2_ymax)
#plot smooth frames
p99 <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_n1, y = ~n1_trend,
data = smooth_frame_n1,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n1,
'</br> Log Copies: ', round(n1_trend, digits = 2),
'</br> Target: N1'),
line = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_lines(x = ~date_vec_n2, y = ~n2_trend,
data = smooth_frame_n2,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n2,
'</br> Log Copies: ', round(n2_trend, digits = 2),
'</br> Target: N2'),
line = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_ribbons(x ~date_vec_n1, ymin = ~n1_ymin, ymax = ~n1_ymax,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n1, #leaving in case we want to change
'</br> Max Log Copies: ', round(n1_ymax, digits = 2),
'</br> Min Log Copies: ', round(n1_ymin, digits = 2),
'</br> Target: N1'),
name = "",
line = list(color = '#1B9E77')) %>%
plotly::add_ribbons(x ~date_vec_n2, ymin = ~n2_ymin, ymax = ~n2_ymax,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_n2, #leaving in case we want to change
'</br> Max Log Copies: ', round(n2_ymax, digits = 2),
'</br> Min Log Copies: ', round(n2_ymin, digits = 2),
'</br> Target: N2'),
name = "",
line = list(color = '#D95F02')) %>%
layout(yaxis = list(title = "Total Log SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
plotly::add_segments(x = as.Date("2020-06-24"),
xend = as.Date("2020-06-24"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "Bars Repoen",
hoverinfo = "text",
text = "Bars Reopen",
showlegend = FALSE,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-07-09"),
xend = as.Date("2020-07-09"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "Mask Mandate",
hoverinfo = "text",
text = "Mask Mandate",
showlegend = FALSE,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-08-20"),
xend = as.Date("2020-08-20"),
y = ~min(n1_ymin), yend = ~max(n1_ymax),
opacity = 0.35,
name = "Classes Begin",
hoverinfo = "text",
text = "Classes Begin",
showlegend = FALSE,
line = list(color = "black", dash = "dash"))
p99smooth_extract <-
plotly::subplot(p99,p1, # plots to combine, top to bottom
nrows = 2,
heights = c(.6,.4), # relative heights of the two plots
shareX = TRUE, # plots will share an X axis
titleY = TRUE
) %>%
# create a vertical "spike line" to compare data across 2 plots
plotly::layout(
xaxis = list(
spikethickness = 1,
spikedash = "dot",
spikecolor = "black",
spikemode = "across+marker",
spikesnap = "cursor"
),
yaxis = list(spikethickness = 0)
)## Warning: Ignoring 1 observationssmooth_extractsave(smooth_extract, file = "./smooth_extract.rda")p1<- p1 %>% plotly::add_segments(x = as.Date("2020-06-24"),
xend = as.Date("2020-06-24"),
y = 0, yend = x_top,
opacity = 0.35,
name = "Bars Repoen",
hoverinfo = "text",
text = "Bars Reopen",
showlegend = TRUE,
data = only_background,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-07-09"),
xend = as.Date("2020-07-09"),
y = 0, yend = x_top,
opacity = 0.35,
name = "Mask Mandate",
hoverinfo = "text",
text = "Mask Mandate",
showlegend = FALSE,
data = only_background,
line = list(color = "black", dash = "dash")) %>%
plotly::add_segments(x = as.Date("2020-08-20"),
xend = as.Date("2020-08-20"),
y = 0, yend = x_top,
opacity = 0.35,
name = "Classes Begin",
hoverinfo = "text",
text = "Classes Begin",
showlegend = FALSE,
data = only_background,
line = list(color = "black", dash = "dash"))